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Abstract 

Avalanche dynamics is an indispensable feature of complex systems. Here we study the 
self-organized critical dynamics of avalanches on scale-free networks with degree exponent 
y through the Bak-Tang-Wiesenfeld (BTW) sandpile model. The threshold height of a node 

1 — Tl 

i is set as k i with < T| < 1, where kj is the degree of node i. Using the branching process 
approach, we obtain the avalanche size and the duration distribution of sand toppling, which 
follow power-laws with exponents x and 8, respectively. They are given as x = (y— 2r\)/ (y— 
1 — r|) and 8 = (y — 1 — r|)/(y— 2) for y < 3 - r|, 3/2 and 2 for y > 3 — r\, respectively. The 
power-law distributions are modified by a logarithmic correction at y = 3 — r\. 
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1 Introduction 

Frequently, complex systems in nature as well as in human society suffer massive 
catastrophes triggered from only a small fraction of their constituents. Unexpected 
epidemic spread of diseases and the power outage in the eastern US of the last 
year are the examples of such avalanche phenomena. Such a cascading dynam- 
ics is not always harmful to us. The information cascades making popular hits of 
books, movies, and albums are good to writers, actors, and singers, respectively. 
Thus it is interesting to understand and predict how those cascades propagate in 
complex system. Recently, the network approach, by which a system is viewed as 
a network consisting of nodes representing its constituents and links interactions 
between them, simplifies complicated details of complex systems. Such a simplifi- 
cation unveils a hidden order such as scale-free behavior in the degree distribution. 
Here degree is the number of links connected to a certain node. The Internet at 
the autonomous system level, the World-Wide Web, social acquaintance networks, 
biological networks, and other many complex networks exhibit power-law degree 
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distributions, pd(k) ~ k~ J . The networks following such power-law degree distri- 
butions are called scale-free (SF) networks [1], where non-negligible fractions of 
hubs, the nodes with extraneously large degrees, exist. 

In this paper, we investigate the avalanche dynamics on such SF networks through 
the Bak-Tang-Wiesenfeld (BTW) sandpile model [2], a prototypical model ex- 
hibiting self-organized criticality (SOC). The study of sandpile dynamics has been 
carried out mostly on regular lattices in the Euclidean space. In the stationary state, 
which can be reached without tuning a parameter, the system exhibit scale-invariant 
features in the power-law form of the avalanche size distribution p a (s) and the du- 
ration or lifetime distribution £(t) as 

p a (s)~s- z and £(t)~r\ (1) 

Recently, Bonabeau has studied the sandpile dynamics on the Erdos-Renyi (ER) 

random networks [3] and found that the avalanche size distribution follows a power 

law with the exponent x ~ 1 .5, consistent with the mean-field solution [4] . Recently, 

Lise and Paczuski [5] studied the Olami-Feder-Christensen model [6] on regular 

ER networks, where degree of each node is uniform but connections are random. 

They found the exponent to be x 1 .65. However, when degree of each node is not 

uniform, they found no criticality in the avalanche size distribution. Note that they 

assumed that the threshold of each node is uniform, whereas degree is not. Here we 

study the BTW sandpile model on SF networks, where the threshold zi of the node 
In 

i is given as k t 1 with ki the degree of i and < T] < 1 . We find that the exponents 
for the avalanche size and the duration distribution depend on the degree exponent 
yas x = (y— 2t|)/(y— 1 —r\) and 8 = (y- 1 -r|)/(y-2) for y < 3 - r\ while, for 
y > 3 — T), they show the same behaviors as the conventional mean-field solutions 
as observed for the ER random networks. 



2 Sandpile model 

We present the dynamic rule of the BTW sandpile model on a given network. 

(1) Each node i is given a prescribed threshold n (< ki). The smallest integer not 
smaller than zi is denoted as \zf\ (\zi] < h). 

(2) At each time step, a grain is added at a randomly chosen node i. The integer- 
valued height of the node i, h t , increases by 1. 

(3) If the height at the node i reaches or exceeds zu then it becomes unstable and 
the \zi] grains at the node topple to its \z[\ randomly chosen adjacent nodes 
among ki ones; 

hi — > hi — \zi] , and hj = hj + 1 for all nodes j which are chosen. 

(4) If this toppling causes any of the adjacent nodes receiving grains to be unsta- 
ble, subsequent topplings follow on those nodes in parallel until there is no 
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unstable node left. This process defines an avalanche. 
(5) Repeat (2)-(4). 

Here the threshold Zi of node i is given as 

a = i, H (0<ii<l), (2) 

which is a generalization of zt = h previously investigated in Ref. [7]. We concen- 
trate on the distributions of (i) the avalanche area A, i.e., the number of distinct 
nodes involving in a given avalanche, (ii) the avalanche size S, i.e., the number of 
toppling events in a given avalanche, and (iii) the duration T of a given avalanche. 



3 Branching process approach 

The mapping of each avalanche to a tree provides a useful way of understanding 
the statistics of avalanche dynamics analytically. For each avalanche event, one can 
draw a corresponding tree: The node where an avalanche is triggered corresponds 
to the originator of the tree and the following nodes to descendants. In the tree struc- 
ture, a descendant born at time t is located away from the originator by distance t 
along the shortest pathway. The tree stops to grow when no further avalanche pro- 
ceeds. Then the ensemble of avalanches can be identified with that of trees grown 
through the branching process. In this mapping, the avalanche duration T is equal 
to the lifetime of the tree minus one, and the avalanche size S differs from the tree 
size only by the number of boundary nodes of the tree, which is relatively small 
when the overall tree size is very large. If one assumes that branching events at dif- 
ferent nodes occur independently and that there is no loop in the tree, the tree size 
and lifetime distribution can be obtained analytically [8,9]. Those distributions are 
expected to share the same asymptotic behaviors with the avalanche size and du- 
ration distribution, respectively, due to the near-equivalence between an avalanche 
and its corresponding tree in their scales as mentioned above. 

In the branching process describing an avalanche, after initial branching into k de- 
scendants with probability qo(k), successive branchings are assumed to occur inde- 
pendently with probability q(k). qo(k) and q(k) may be different in general, but the 
statistics of the overall size and duration of an avalanche is determined dominantly 
by q{k) . We checked also numerically the case where a new grain is added to a node 
with the probability proportional to the degree of that node, which gives different 
qo(k) from that in the case where a new grain is added randomly, and found that the 
nature of the avalanche dynamics is the same in both cases. Thus, for simplicity, 
we consider the branching process where every branching occur with probability 
q(k). For the BTW model in the Euclidean space, where the threshold zi of node i is 
equal to its degree ku q(k) has a finite cut-off such that q(k) = for k > Zi = const, 
because the degree of each node is uniform and finite. Consequently, the exponents 
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of the avalanche size and the duration distributions in Eq. (1) come out to be the 
so-called mean-field values; x = 3/2 and 8 = 2 [8,9]. These results are known to 
hold for the BTW model on regular lattices with dimensions larger than 4 [4]. Note 
that when dimension is smaller than 4, the branching process approach cannot be 
applied, so that the values of the exponents x and 8 would not be trivial. 

In SF networks, avalanches usually do not form loops, generating tree-structures: 
According to the numerical simulations of the BTW model for the case of zi = h 
on SF networks [7], the statistics of the two quantities A and S are nearly equal 
when they are large: For example, the maximum area and size (A max , S mRX ) among 
avalanches are (5127, 5128), (12058, 12059) and (19692, 19692) for scale-free 
networks with y= 2.01, 3.0, and oo, respectively. The fact that A and S are almost 
the same implies that the avalanche structure can be treated as a tree. From now on, 
we shall not distinguish A and 5, denoted by s. Thus it is valid to use the branching 
process approach to understand the avalanche dynamics on SF networks. 

We study the BTW model on SF networks with the degree exponent y and the 
threshold given as Eq. (2). The branching probability q(k) consists of two factors, 
that is, q(k) = q\(k)q2(k), where qi(k) is the probability that the threshold zi of 
node i is in the range k — 1 < Zi <k and qi{k) is the probability that the total 
number of grains at the node reaches or exceeds the threshold. If Zi = f{h) with 
f{x) a monotonic increasing function of x satisfying f(x) < x for all x > 1, the 
condition of k — 1 < Zi < k implies that q\{k) is nothing but the probability that 
a node i connected to the one end of a randomly chosen edge has its degree ki 

in the region {f~\k- l),f~ l (k)], and thus qi (k) = I,lQ\- l)l+ i *Pd(*)/ '<*>, 
where |*J is the largest integer not larger than x. Notice that T^ = iqi(k) = 1 and 
qi(k) ~ A:( 1 -Y+ T l)/( 1 - r l) for large k if f(x) -x 1 ^ (0 < T| < 1) for large x. q 2 (k) is 
the probability that the node i has height k — 1 at the moment of receiving a grain 
from one of its neighbors. We have checked numerically that a typical height of 
node is absent, so that all possible k values 0, 1, . . . ,k — 1 are equally likely [10]. 
Thus we set q2{k) = l/k. As a result, the branching probability q{k) for large k is 
given asymptotically as 

*(*) = i*i (*)-*-* {i = T^)- ( 3 ) 

When zi = h or r\ = 0, y 1 is reduced to y. Since we are interested in the case of y > 2 
and < r| < 1, y > 0. 

Using the independence of the branchings from different parent-nodes, one can de- 
rive the following self-consistent relation for the tree size distribution p(s) as [8,9] 



pCs) = <7(0)8 M £ £ ••• Y,p(si)p(s2)...p(s k )%* =iSiiS _ v (4) 

k=l .si = 1a'2 = 1 s^ = \ 
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This relation can be written in a more compact form by introducing the generating 
functions, T{y) = l^i pW and Q,(to) = 1^=0^)^ as 

<P{y)=yd{<P{y)). (5) 

Then to = •£ (y) is obtained by inverting y = <E ~ 1 (to) = to/Q, (to) . 

The average size (s) of a finite tree can be obtained easily from the generating 
functions. 

finite 

( s )= = *'(1), (6) 

.5=1 

where fp'(y) = d<P{y)/dy. Using the relation, Eq. (5), we obtain 

{) [) i-o/(*(i))' (j 

where again Q,'(to) = dQ,((Q)/d(Q. 

The distribution of duration, z'.e., the lifetime of the tree can be evaluated simi- 
larly [8,9]. Let r{t) be the probability that a branching process stops at or prior to 
time t. Then following the similar steps leading to Eq. (4), i.e., r(t) = Er=o^D"( ? — 
l)] k , one has 

r(t) = d(r(t-l)). (8) 
For large t, r(t) comes close to 1. One can obtain (0 = r(t — 1) by solving doo/dt ~ 
r{t) —r(t — \) = <2(to) — (0. Then the lifetime distribution £(t) is obtained through 
£(t) = r(t) - r(t - 1) ~ rfto/J?. 



4 Avalanche size and duration distribution 



The growth of a tree depends on the average number of branches defined as 

oo 

C=Y,kq(k). (9) 

k=\ 

When C > 1 (C < 1), a tree can (cannot) grow infinitely in a probabilistic sense. 
Thus the case of C = 1 is a critical point for the growth of a tree. One can see 
that for any branching process with q(k) = (\/k)q\(k) (k > 1) and YX=\ = 1' 
the average number of branches C is always 1, independent of detailed structural 
properties of networks. Therefore our assumption q2{k) — l/k corresponds to the 
condition for the self-organized criticality (SOC) of the sandpile model. 

The inverse function !P _1 ((o) satisfies !P _1 (1) = 1. When C = 1, the first-order 
derivative 32>~ 1 (to) /da> at (0 = 1 is zero and thus <B (y) becomes singular at y = 1. 
2>(v) is expanded around y = 1 as <p{y) ~ 1 — b{\ —y)^ with constant b and < (])< 
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1 . Then the asymptotic behavior of the avalanche size distribution p(s) for large s is 
given by p(s) ~ s^ -1 , because if a series Y,7=o a sy s with the radius of convergence 
1 has the asymptotic behavior 

oo 

£ a s f ~ (1 -y^ as y -> 1, then a s . ~ as 5 -> oo. (10) 

5=0 



The functional form of the branching probability g(fc) determines the singularity of 
fp (y). To illustrate this, we first consider a simple case that 



q(k) 



I -a (k = 0), 

a (A: = 2), (11) 

(otherwise), 



where < a < 1. Then the average number of branches C = Y,kq(k) = 2a and 
the generating function Q (a)) = Lr=o l(k)d> k = \ —a + a(D 2 . Using the relations of 
y = cq/q,(cq) and a> = 4> (y), it is obtained that 

The value of 4>(1) = X^fpOs') is given as 

, 1N l-|l-2a| f 1 for 0<a<± (C < 1), 

which means that when 1/2 < a < 1 (C > 1), a tree can grow infinitely with prob- 
ability 1 — fp(l) = (2a — l)/a, and the critical point is located at a c = 1/2. Near 




y = 1, <P(y) « 1 — a/2(1 — y) from Eq. (12), leading to (]) = 1/2. Then, the avalanche 
size distribution p(s) behaves as p(s) ~ s~ 3//2 . On the other hand, using Eq. (7), 

(a < a c ), 



(s) = <p'{\) = { 2{ac l a) (14) 
4^) {a>a L ). 



Even for the case that q{k) has a finite cut-off larger than 2 or decays exponen- 
tially, the above result holds. This is the conventional mean-field solution for the 
avalanche size distribution [4,8,9] and has been shown to hold for the BTW model 
on the ER random networks [3]. 

When q(k) decays slowly as in Eq. (3), however, its generating function d((D) is 
singular at to = 1. For q(k) in Eq. (3), the expansion of d(K>) around to = 1 is given 
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as 



q,(q>) ~ 1 — (1 — co) + 



Ai(l-a))^- 1 (2<y<y c ), 
-A 2 (l-(o) 2 ln(l-(o) (y=y c ), (15) 
A 3 (l-(o) 2 (Y>Yc), 



where A,'s are constants, y 7 is given in Eq. (3), and y c = 3 — T|. The derivation of 
the logarithmic correction for the case of y = y c can be found in [11]. Note that the 
singular term (1 — O))"^ -1 is the second leading term of 1 — Q,((xi) for y < y c . Using 
the relation fP _1 (co) = gd/o(g)) in Eq. (5), the behavior of <p{y) around y = 1 is 
obtained for each region of y from Eq. (15), and in turn, using Eq. (10), p(s) for 
s — > oo. We find that 



p(s) ~ < 



s _( y -2t 1 )/(y-1-t 1 ) ( 2 <y< Yc ) ) 
s-^ 2 {\ns)- 1 / 2 (y=y c ), 

S- 3 /2 (Y>Yc)- 



(16) 



Thus, the exponent! is given asx = (y-2r\) / (y- 1 — rj) for 2 < y< y c and x = 3/2 
fory> Yc- 

Also obtained is r{t) from Eq. (15) by using Eq. (8). The duration distribution £(t), 
which is the derivative of r(t), is found to be 



r (T-i-ri)/(Y-2) (2< Y <Yc), 
t-^lnt)^ ( Y =Yc), 
^ 2 (Y>Yc)- 



(17) 



That is, the exponent 8 is given as 8 = (y — 1 — T|)/(y — 2) for 2 < y < y c and 8 = 2 
fory> y c . 



5 Conclusion 



We have studied the BTW sandpile model on SF networks with the degree exponent 
yto understand the avalanche dynamics in complex systems. The main results are 
the avalanche size and duration distribution. The exponents x and 8 increase with 
increasing y, implying that the hubs play a role of reservoir, that is, sustain large 
amount of grains to make the SF network resilient under avalanche dynamics. This 
is reminiscent of the structural resilience of the SF network under random removal 
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of nodes for y < 3 [12,13,14]. We also checked the case where the threshold Zi 
contains noise in the way that n = C,iki with being distributed uniformly in [0,1]. 
We find that such a variation does not change the nature of the avalanche dynamics. 
However, when the threshold is given in terms of a quantity other than degree, 
e.g., load, the corresponding avalanche dynamics has no reason to follow the same 
statistics as studied in this paper, which remains further works. 
This work is supported by the KOSEF Grant No. R 14-2002-059-0 1000-0 in the 
ABRL program. 
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